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Abstract 

We  develop  a  general  class  of  approximations  of  mean-spherical  (MSA)  type  as  a 
method  for  studying  lattice  percolation  problems.  We  review  the  MSA  and  test  certain 
extensions  of  it  on  lattice  spin  models.  The  relations  between  thermal  spin  models  and 
percolation  models  are  then  reviewed  in  order  to  identify  natural  extensions  of  the  MSA 
to  percolation  models.  These  extensions  are  used  to  treat  both  site  and  bond  percolation 
models.  In  one  ’low-density’  formulation  of  MSA,  the  threshold  for  bond  percolation  on  a 
lattice  is  found  to  equal  the  value  at  the  origin  of  the  corresponding  lattice  Green’s  function. 
This  result  is  extremely  accurate  for  all  lattices  studied,  and  in  all  space  dimensions  d  £ 
3.  An  accurate  treatment  is  also  given  of  the  general  site-bond  problem.  The  entire 
percolation  locus  for  this  problem  agrees  very  closely  with  the  results  of  simulation.  We 
also  introduce  a  new  method  for  studying  percolation  transitions  which  is  a  hybrid  of 
the  Kikuchi  cluster  approximation  scheme  and  the  MSA.  The  method  is  shown  to  give 
extremely  good  values  for  percolation  thresholds  while  preserving  the  valuable  features  of 
the  standard  MSA.  In  particular,  it  provides  a  convenient  means  of  computing  the  pair 
connectedness  function.  We  also  explore  extensions  of  our  approximations  to  treat  directed 


site  and  bond  percolation. 


1.  Introduction 


Percolation  is  a  phenomenon  defined  by  the  formation  of  macroscopic  clusters  in  a 
many*body  system,  given  a  criterion  for  pairwise  connectedness.  It  has  been  studied  in 
recent  years  from  at  least  two  different  points  of  view.  Those  studying  percolation  as  a 
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particularly  accessible,  geometrical  analog  of  a  phase  transition  seek  an  accurate  method 
to  locate  the  percolation  threshold  and  to  study  the  scaling  behavior  of  physical  quantities 
in  its  vicinity.  On  the  other  hand,  those  who  want  to  study  the  influence  of  percolation 
on  bulk  material  properties  and  transport  processes  in  disordered  materials  and  liquids 
want  an  account  of  various  key  quantities  that  describe  clustering,  such  as  the  average 
coordination  number  of  a  particle  and  the  mean  cluster  Bize,  throughout  the  entire  range 
of  system  parameters.  An  approach  that  addresses  both  classes  of  problems  is  that  based  on 
the  solution  of  integral  equation  approximations  for  the  two-point  connectedness  function. 
Several  such  integral  equations  have  been  derived  and  solved,  both  for  continuum  and 
lattice  percolation.1-8 

One  such  integral  equation  approximation  which  has  been  studied  for  various  con¬ 
tinuum  percolation  models  »  the  Ornstein-Zernike  equation  with  closure  provided  by  the 
mean  spherical  approximation  (MSA).  This  approximation  was  originally  developed9  for 
nearest-neighbor  thermal  6pin  models  and  their  equivalent  lattice-gas  analogs.  It  was  sub¬ 
sequently  generalized10  to  longer-range  lattice  models  and  to  continuum  fluid  models.  The 
MSA  for  a  continuum  system  of  particles  at  thermal  equilibrium  can  be  defined  by  the 
equations 


M*12)  —  —  1  *12  <  a 


«(*12)  =  -0v(xi2)  *12  >  O 


(1.1) 


Here,  h(xi2)  and  c(*j2)  are,  respectively,  the  full  correlation  function  and  the  direct  corre¬ 
lation  function  between  particles  at  the  vector  positions  *1  and  x2.  Here,  *12  =  |*i  —  *2! 
is  the  scalar  distance  between  such  particles,  17(3:12)  is  their  pairwise  potential  energy,  and 
/?  =  (kT)-1,  where  T  is  the  absolute  temperature,  and  k  is  Boltzmann’s  constant.  The 
distance  a  is  the  distance  of  closest  approach  fixed  by  a  hard-core  potential  which  is  part 


□ 

□ 


of  the  pairwise  interaction.  The  MSA  has  been  used  profitably11  to  treat  a  wide  variety  of 


thermal  problems.  A  major  benefit  of  this  approximation  is  that  the  resulting  equations, 

C.Qlli 

e.g.,  for  the  equation  of  state,  can  be  solved  analytically  for  manv  models  of  interest. 
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Recently4  the  MSA  was  extended  to  relate  the  two-point  connectedness  function  hc(x ) 
and  the  direct  connectedness  function  cc(x)  of  the  random  sphere  percolation  model.  This 
is  a  model  of  randomly  placed  spheres  which  are  taken  to  be  directly  connected  if  their 
pairwise  separation  is  less  than  a  fixed  distance  a,  the  sphere  diameter.  For  this  model, 
the  MSA  is  equivalent  to  the  Percus-Yevick,  which  was  first  suggested  by  Coniglio  et  als 
in  work  that  extended  the  cluster-  expansion  treatment  of  he(x )  by  Hill.6  More  general 
models  have  been  introduced  by  considering  systems  of  particles  in  thermal  equilibrium 
interacting  via  a  pair  potential  <j>(x),  which  introduces  correlation  between  the  centers 
of  particles  (if  one  sets  <t>{x)  —  0  the  model  reduces  to  the  random  sphere  percolation 
model).  A  number  of  such  continuum  systems  of  interest  can  be  solved  exactly  in  the 
MSA.T“®  *****  It  is  also  natural  to  apply  the  MSA  to  the  functions  hc{x)  and  cc(*)  in 
lattice  percolation  models.  A  recent  study  of  MSA  for  lattice  Bite  percolation  by  Hoye  and 
Stell6  found  that  both  the  percolation  threshold  and  critical  exponents  for  percolation  on 
certain  lattices  in  three  and  higher  dimensions  are  predicted  accurately.  •  In  this  paper, 
we  extend  that  study  by  using  certain  general  methods  of  MSA  type  to  investigate  lattice 
bond  and  site-bond  percolation.  The  critical  exponents  are  found  to  be  the  same  as  in 
site  percolation,  as  one  expects,  while  the  bond  percolation  threshold  is  predicted  with 
remarkable  accuracy  on  all  lattices  studied,  and  in  all  dimensions  d  >  3.  In  this  paper, 
we  will  work  with  a  very  general  class  of  approximations  of  mean  spherical  type.  These 
approximations,  in  general,  define  a  closure  of  the  Ornstein-Zernike  equation  by  providing 
a  pair  of  assumptions  corresponding  to  the  two  given  in  (1.1).  In  general: 

1.  The  volume  surrounding  each  lattice  site  is  divided  in  two  by  choosing  a  sphere  that 
surrounds  that  site.  The  value  of  the  pair  correlation  function  (either  h(x)  or  c(*)  can  be 
used),  giving  the  interaction  of  a  site  with  other  sites  inside  the  chosen  sphere,  is  provided 
explicitly,  either  L,  the  constraints  of  the  model  itself,  or  by  some  other  approximation. 
Many  approximation  schemes  provide  such  values  for  correlation  functions  at  small  sepa¬ 
ration. 
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2.  The  interaction  with  lattice  sites  outside  the  sphere  is  described  by  giving  an  approx¬ 
imate  form,  for  large  separation,  of  the  direct  correlation  function  c(x).  In  the  original 
MSA,  as  given  by  (1.1),  the  radius  of  the  sphere  in  the  scheme  just  described  is  taken  to 
be  the  exclusion  radius,  defined  as  the  range  of  the  hard-core  part  of  the  potential,  and 
the  closure  is  specified  at  small  separation  by  giving  the  values  of  h(se)  inside  the  exclusion 
sphere.  We  note  that  the  scheme  just  outlined  has  obvious  similarities  with  the  effective 
medium  approximation  (EMA). 15,11  It  also  has  basic  differences;  for  example,  it  gives,  in 
general,  critical  exponents  that  are  non-classical. 

This  paper  is  organized  as  follows:  in  section  2,  we  review  the  general  procedure  for 
solving  the  MSA  for  thermal  lattice  models.  We  illustrate  the  use  of  MSA-like  approxima¬ 
tions  by  applying  these  to  the  Ising  model.  In  section  3,  we  discuss  the  relation  between 
thermal  models  and  percolation  models  in  order  to  identify  natural  extensions  of  the  MSA 
framework  to  percolation.  We  review  the  earlier  treatment  of  site  percolation.  Appropriate 
variants  of  the  MSA  are  then  applied  to  both  bond  and  site-bond  percolation.  In  Section 
4,  we  discuss  the  use  of  the  Kikuchi  cluster  approximation  to  obtain  the  short-distance 
values  of  h(x)  needed  to  complete  an  MSA-like  closure.  In  Section  5,  we  apply  MSA- 
like  approximations  to  directed  site  and  bond  percolation  models.  Section  6  discusses  the 
present  limitations  of  the  approximations  discussed  in  this  paper,  and  gives  suggestions 
for  further  research.  Section  7  summarizes  our  conclusions. 

2.  General  MSA  Formalism 

In  this  section,  we  review  the  general  Ornstein-Zernike  formalism  for  lattice  models. 
This  formalism  is  identical  for  thermal  and  percolation  problems.  We  illustrate  its  use  by 
applying  it  to  the  nearest-neighbor  Ising  model,  as  described  in  lattice-gas  terminology. 

The  Ornstein-Zernike  equation  is 

M*«)  =  ®(*12 )  +  P^j  c(xis)Mx»2)  (2.1) 

*3 

In  this  section  we  use  the  terminology  of  thermal  physics,  in  which  h(zi?)  and  c(*i2)  are, 
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respectively,  the  full  pair  correlation  function  and  direct  pair  correlation  function.  The 
sum  in  2.1  is  over  all  lattice  sites.  We  note  that  the  function  h(x)  and  other  functions 
defined  on  a  lattice  depend  upon  a  vector  which  we  denote  as  x  or  x,j ,  this  is  the  separation 
of  sites  i  and  j.  An  equation  identical  to  (2.1),  but  with  different  boundary  conditions, 
governs  the  relation  between  the  pair  connectedness  function,  given  by  <7c(x)  =  hc(x)  +  1, 
and  direct  connectedness  function  cc(x),  in  a  percolation  model.  This  will  be  discussed  in 
detail  in  the  next  section. 

We  indicate  the  Fourier  transform  of  a  function,  e.g.  c(x),  by  placing  a  caret  over  the 
corresponding  function  symbol.  Thus  we  have 

c(*)=^]J  Jnk)e-"  'd'k  (2.2) 

Here  the  integral  over  wave  number  k  is  over  a  single  Brillouin  zone  of  the  lattice.  Also,  ft  is 
the  volume  of  the  Wigner-Seitz  cell  associated  to  a  single  lattice  site;  it  is  here  normalized 
to  unity.  Taking  the  Fourier  transform  of  both  sides  of  (2.1)  allows  an  algebraic  solution 
for  h(k)  in  terms  of  c(k): 

w 

1  —  pc(k) 

The  strategy  for  solving  (2.1)  is  this:  since  c(x)  is  assumed  to  be  short-ranged,  c(k)  can  be 
written  explicitly  in  terms  of  a  small  number  of  unknown  constants,  namely  the  different 
values  of  c(x)  which  are  non-negligible.  Substituting  an  assumed  form  for  c(fc)  in  (2.3)  and 
Fourier-  transforming  then  gives  a  solution  for  h(x).  Boundary  conditions  and  constraint 
conditions  specify  the  remaining  constants  in  the  solution.  These  conditions  contain  the 
physical  assumptions  specific  to  the  model  being  studied. 

Specifically,  for  the  problems  of  interest  in  this  section,  we  assume  that  c(x12)  is 

nonzero  only  if  x12  is  either  zero  or  equal  to  a  nearest-neighbor  lattice  vector,  in  which 

cases  c(x)  takes  the  values  Co  and  c2,  respectively.  Substituting  this  information  in  the 

transform  inverse  to  (2.2)  then  gives,  for  hypercubic  lattices 

d 

e(k)  =  co  +  2cj  ^  cos  k{a  (2.4) 

«=i 


m 
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Here  a  is  the  lattice  spacing  and  d  is  the  dimensionality  of  the  lattice.For  the  sake  of 
algebraic  simplicity,  we  specialize  our  discussion  to  hypercubic  lattices  and  will  work,  unless 
otherwise  stated,  with  the  three-dimensional  simple  cubic  (SC)  lattice.  The  discussion  of 
thi6  paper  applies,  however,  to  general  Bravais  lattices. 

Taking  the  Fourier  transform  of  (2.3)  and  adding  to  both  sides  the  identity 


J  r: 


“pc(fe)  .-«*• 

pc{k) 


xd3k 


(2.5) 


then  gives 

=  (20) 

Using  (2.4)  for  c(k)  allows  us  to  rewrite  the  denominator  of  the  RHS  as 


1  —  pc(fc)  =  zpc  i 


d 

(1  -I-  K 2)  —  d~l  cos  kja 

»=i 


(2.7) 


Here  z  is  the  coordination  number  for  the  lattice,  and  the  quantity  K2  is  given  by  the 
expression 

K 2  =  zk 2 

_  1  -  pco  -  ZPC1  (2-8) 

pci 

Substituting  (2.7)  into  (2.6)  gives 


£*,o  +  pM*) 


J _ f  «~a- 

zpc\  (2tt)3  J  (1  +  k2  )  —  d~l  cos  kia 
zG(x) 
zpci 


(2.9) 


This  last  equation  defines  G(z),  which  can  be  identified  as  the  Green’s  function  of  the  lattice 
version  of  the  Helmholtz  wave  equation  with  wavenumber  K 2  =  zk2.  This  function  occurs 
frequently  in  mathematical  physics  and  has  been  tabulated.16’10  We  note  in  particular  that 
for  k 2  =  0,  the  function  zG(x)  is  the  generating  function  for  random  walks  on  the  lattice 
being  studied. 
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We  digress  briefly  to  justify  our  identification  of  k,  as  defined  by  (2.8),  with  the  inverse 
correlation  length.  For  convenience,  define  x(*)  =  5*,o  +  ph{z).  Then  we  have 


X(0)  _  l-fM(Q) 
X(k)  1  +  ph(k) 
1  —  pc(k) 
~  1  -  pc( 0) 


(2.10) 


where  the  second  step  uses  the  Ornstein-Zernike  equation  (2.3).  Expanding  the  RHS  of 
(2.10)  in  powers  of  k 2  and  substituting  (2.4)  gives 


1  4-  A2  (fca)2  +  0[(fca)4] 


with  k2  the  norm  squared  of  the  vector  k  and  A2  defined  by 


(2.11) 


A2  = 


1  -  pc(0) 


(2.12) 


For  simple  hypercubic  lattices,  A2  can  equivalently  be  identified  as17’18 

A  E*  ***(*) 

2"  £.x(*) 


(2.13) 


That  is,  A2  is  the  second  spherical  moment  of  the  quantity  x(*)  divided  by  its  zeroth 
moment.  To  see  this,  write  x(*)  in  terms  of  its  Fourier  transform,  replace  the  factors  of  * 
with  ib-derivatives,  and  note  that  (2.11)  is  a  Taylor  series.  Direct  comparison  then  shows 
that  zA2  =  k“2,  with  k 2  defined  by  (2.8).  The  singular  part  of  the  susceptibility  near  a 
critical  point  is  proportional  to  the  volume  integral  of  h(z),  i.e.,  to  h(0).  According  to 
(2.3),  we  can  then  identify  such  a  point  by  the  condition 


pc(0)  =  1 


(2.14) 


Note  that  by  (2.8),  this  is  equivalent  to  the  condition  k 2  =  0.  Also,  it  is  reasonable  to 
identify  A2  with  the  square  of  the  correlation  length  and  thus  identify  k~2  =  £2.  The 
correlation  length  defined  in  this  manner  has  the  same  critical  behavior  as  that  based 
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on  the  assumption  that  h(x)  decays  exponentially;  the  former  definition  is  however  more 
generally  applicable.  For  example,  even  when  h(x)  has  algebraic  decay,  as  it  will  when  c(z) 
does,  this  definition  is  still  applicable,  providing  that  the  second  moment  of  c(x)  exists.11*1 
Several  further  relations  are  useful  for  the  general  discussion  of  critical  behavior  as 
described  by  the  Ornstein-Zemike  equation.  To  develop  these,  we  write  Go  and  G\ ,  respec¬ 
tively,  for  the  values  of  the  Green’s  function  G(x),  defined  by  (2.9),  evaluated  for  a  lattice 
displacement  vector  equal,  respectively,  to  zero,  and  to  the  vector  difference  between  a  pair 
of  nearest  neighbors.  We  can  specify  these  displacement  vectors  as  x  =  0,  and  |x|  =  a, 
respectively.  Note  that  this  causes  no  ambiguity  because,  by  symmetry,  G(x)  takes  the 
same  value  for  all  nearest  neighbor  displacement  vectors.  In  the  special  case  that  k 2  =  0, 
which  defines  the  critical  point,  we  write  these  same  quantities  as  G0  and  G\ ,  respectively. 
Similarly,  the  values  of  h(x)  at  x  =  0  and  |z|  =  a  are  written  ho  and  hi,  respectively. 
Setting  |z|  =  a  in  (2.9)  then  gives 

—  -  phi  ■  (2.15) 

zpci 

Dividing  (2.9)  for  x  =  |aj  by  the  same  equation,  with  *  =  0,  gives 

zGi  Phi 

zGo  1  +  pho  K  ' 

Finally,  we  can  relate  Go  and  G i,  for  a  SC  lattice,  by  using  the  identity 

n  f - m  zGt  (2.,7) 

(2tt)s  7(14-  k2)  -  d~l  53*=i  cos  k*a 

(A  similar  identity  holds  for  any  lattice.)  This  follows  directly  by  using  the  symmetry  of 
a  Bravais  lattice  and  the  definition  of  Gi.  This  gives  directly,  using  (2.5)  and  (2.9) 


(1  +  k?)zGo  —  zG  i  —  1 


(2.18) 


This  relation  can  be  used  to  eliminate  G\  from  (2.16)  and  give  a  basic  relation  for  the 


critical  point: 


zGp  -  1  phi 

zQo  1  4-  pho 


(2.19) 
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We  note  that  the  lattice-gas  density  p  occurs  explicitly  only  in  the  combinations  pc0, 
pc\.  Two  further  constraints  must  now  be  supplied  to  completely  specify  an  MSA-like 
closure  of  the  Ornstein-Zernike  equation.  In  standard  applications  of  the  MSA  to  pair 
correlation  functions,  these  tie  chosen  to  be  the  vanishing  of  the  correlation  h(x)  inside 
the  interaction  hard  core,  and  the  equation  c(®)  =  —/3v(x)  outside  this  core.  Here  v(z) 
is  the  interparticle  potential  and  /?  =  1/kT.  The  second  of  these  equations  is  a  ’’linear 
response”  high-temperature  approximation;  it  is  exact  to  first  order  in  /3  and  the  potential 
tr(z).  For  the  Ising  model,  this  implies  h0  =  —  1  and  c\  =  K  where  K  =  -(3  Jlg  and  Jlg 
is  the  lattice-gas  coupling  constant,  related  to  the  Ising  model  coupling  constant  J ising 
by 

Jlg  =  4  J ising  (2.20) 

Substituting  these  relations  into  (2.15)  and  (2.19)  then  gives  the  condition  for  criticality 

p(l  -  p)Kcrit  =  Go  (2.21) 

Because  we  consider  the  high  temperature,  zero-field,  Ising  model  only,  we  have,  for  the 
lattice-gas  density,  p  =  0.5.  Thus 


Kcrit  —  4Cm>  ~  1.08 


(2.22) 


for  the  three  dimensional  simple  cubic  (SC)  lattice.  This  should  be  compared  with  the 
value  Kent  =  *840  given  by  the  second  Bethe-Peirls  approximation  for  this  lattice,  and  also 
with  the  value  Kerit  =  .918  given  by  high-temperature  series  analysis,19  which  represents 
the  most  precise  means  available  for  estimating  such  quantities. 

We  derive  a  variant  of  this  approximation  if  we  recall  the  definition  of  c(®)  as  the 
direct  correlation  function  and  write  the  low-density  approximation 


c(x)  =  f(x) 

=  e~K  —  1 


(2.23) 
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This  approximation  is  exact  to  first  order  in  the  lattice-gas  density;  i.e.,  it  is  the  first  term 
in  the  Mayer  expansion.  Using  cj  from  (2.23)  gives  the  criticality  condition 

Kerit  =  —  ln[l  ■+■  40o]  (2.24) 

which  implies  K^t  =  .732.  Thus  for  the  SC  Ising  model  in  three  dimensions,  the  estimates 
for  the  critical  point  given  by  the  high-temperature  approximation  (2.22)  and  the  low- 
density  approximation  (2.24)  are  of  comparable  accuracy.  It  is  (2.22)  (or  in  general  the 
constraint  c(x)  =  —f3v(x))  that  has  come  to  be  called  the  mean-  spherical  approximation 
(because  (2.22)  proves  to  be  exact  in  the  mean-spherical  and  spherical  models  of  a  magnetic 
system.)  As  we  shall  see,  in  some  percolation  problems,  in  which  h(x)  is  replaced  by  the 
pair  connectedness  function  and  c(x)  by  the  direct  connectedness  function,  it  turns  out 
that  (2.23)  (or  in  general  c(x)  =  exp[— fiv(x)]  —  1  )  appears  to  be  the  more  natural  and 
generally  useful  approximation.  This  is  also  found  to  be  the  case  in  studies  of  continuum 
percolation.6 

3.  The  MSA  for  Site  and  Bond  Percolation 
In  this  section,  we  will  apply  the  Ornstein-Zernike  formalism  of  Section  2  to  lattice 
percolation  models.  First,  the  case  of  site  percolation  is  reviewed.  We  show  that  a  nat¬ 
ural  approximation  for  bond  percolation  reproduces  an  analytic  formula  for  the  critical 
point  that  was  already  shown  empirically,  by  Sahimi25  and  co-workers,  to  be  remarkably 
accurate.  Finally,  we  give  sun  approximation  of  MSA  type  for  the  general  case  of  random 
site-bond  percolation  and  reproduce  the  complete  percolation  locus  for  that  model. 

The  formalism  of  Section  2  can  be  directly  applied  to  percolation  models  because  these 
satisfy  an  Ornstein-Zernike  equation  of  the  form  (2.1).  In  particular,  one  has2 

ffc(*12)  =  Cc(*12)  +  Ce(xj j )pc(® J2 )  (3.1) 

*3 

■with  <7c(*i2)  and  ce(*i2)»  respectively,  being  the  connectedness  function  and  direct  con¬ 
nectedness  function,  respectively.  The  Ornstein-Zernike  equation  can  be  taken  to  be  a 
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definition  of  e(x),  and  thus  has  no  content,  per  se.  However,  we  focus  here  on  the  ’deriva¬ 
tion’  of  this  equation  from  other  formulations  that  permit  a  more  direct  estimate  of  the 
value  of  c(x)  and  h(x)  for  small  separation.  This  is  valuable  because,  in  any  case,  one 
needs  an  independent  representation  for  cc(x ),  or  an  independent  second  relation  between 
he(x)  and  ce(x)  in  order  to  have  an  closed  set  of  equations  for  these  quantities.  Below,  we 
shall  refer  to  both  the  p  and  (3  expansions,  which  facilitate  estimation  of  cc(x)  for  small 
argument.  In  the  case  of  random  or  uncorrelated  percolation,  hc  and  cc  are  temperature- 
independent,  so  that  only  the  p  expansion  is  available.  We  focus  in  particular  on  mappings 
of  percolation  models  onto  limiting  cases  of  thermal  models;  these  allow  us  to  draw  on  our 
experience  with  MSA-like  approximations  for  the  latter.  Such  mappings  allow  us  to  exploit 
the  machinery  of  liquid-state  physics.  Also,  they  are  essential  when  thermal  correlations 
are  imposed  between  the  sites  or  bonds  of  a  percolation  model. 

In  the  absence  of  such  correlations,  one  can  calculate  series  for  ec(x)  and  pc(x)  in 
purely  graph-theoretic  terms,  using  a  formulation  due  to  Essam20  in  terms  of  self-avoiding 
walks.  For  example,  for  pure  site  percolation  one  has 

St(x)  =  Y,  d(GK(C|  (3.2) 

G 

where  the  sum  is  over  one-irreducible,  two-rooted  subgraphs  G  of  the  lattice  being  studied, 
v(G)  is  the  number  of  vertices  in  the  graph  G ,  and  d(G )  is  a  purely  combinatoric  quantity 
depending  only  on  the  graph  G.  The  function  cc(x)  is  given  by  a  similar  expression,  but 
with  the  sum  restricted  to  non-nodal  graphs.20 

The  techniques  of  liquid-state  physics,  such  as  integral  equations,  Mayer  expansions, 
etc.,  are  most  naturally  formulated  in  the  continuum,  or  more  precisely,  in  lattice-free 
language  that  applies  to  both  continuum  and  lattice  models.  The  same  is  true  of  the 
formulation  of  percolation  theory  best  suited  to  adapting  these  techniques,  that  due  to 
Hill.1  We  sketch  this  first.  Specifically,  we  write  the  Boltzmann  factor  for  a  thermal  model 
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as  a  sum  of  two  terms: 


exp[-v(x12)]  =  |exp[-^v(®12)]pb(xi2)  J  +  j  exp[-/3v(x12)][l  -p6(x12)] 

=  e+(xa2)  +  e*(x12) 

This  induces  a  corresponding  separation  of  the  Mayer  function 

/(*12)  =  /+(®12)  +  /*(xi2)  (3.4) 

with  /+  =  e+  and  /*  =  e*  —  1.  The  function  p&(xi2)  defines  the  separation-dependent 
probability  of  a  bond  between  two  particles.  Its  choice  is  dictated  by  the  physical  phe¬ 
nomenon  to  be  modelled.  The  first  term  in  (3.2)  is  identified  with  the  particles  being 
directly  connected,  the  second  with  them  not  being  directly  connected.  Substitute  the 
sum  (3.4)  for  each  Mayer  bond  in  the  virial  expansions  of  fi(xJ2)  and  c(xi2),  and  expand 
each  Mayer  graph  into  subgraphs  whose  lines  correspond  to  /+  or  f*  bonds.  Define2,21 
the  connectedness  function  gc(x J2)  to  be  the  sum  of  all  such  subgraphs  in  the  expansion 
of  h(x)  in  which  the  root  points  are  joined  by  a  chain  of  /+ -bonds;  the  blocking  func¬ 
tion  <7&(*i2)  is  the  sum  of  all  the  remaining  subgraphs.  Similarly,  one  defines  the  direct 
connectedness  function  cc(x)  to  be  the  sum  of  the  corresponding  subgraphs  contributing 
to  c(x);  this  is  equivalent  to  the  set  of  subgraphs  contributing  to  </c(x)  that  in  addition 
have  no  nodal  points.  These  definitions  are  compatible  with  (3.1).  Thus  the  correlation 
function  has  been  written  as  the  sum 

M*12)  =  0c(*12)  +  36(*12)  (3-5) 

of  the  two-point  connectedness  function  and  two-p*>int  blocking  function.  One  can  show21 
that  0c(*i2)  as  defined  above  formally  is  in  fact  the  two  point  connectedness  function 
for  a  many-body  system  of  particles  having  correlation  function  h(x)  and  density  p,  and 
being  pairwise-  connected  with  separation-dependent  probability  pt(x).  The  connectedness 
function  £e(x12)  is  the  probability  distribution  associated  with  finding  particles  at  xj  and 
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*a  in  the  same  connected  cluster.  Similarly,  <7fc(*i2),  the  two-point  blocking  function,  is 
the  corresponding  probability  distribution  associated  with  the  particles  at  x\  and  x2  being 
in  different  clusters.  We  note  that  our  probability-density  definition  of  pc(*)  end  the 
corresponding  definition  of  ce(x)  induced  by  (3.1)  are  not  dependent  on  the  density  series 
expansions  of  these  quantities;  the  latter  are  not  fully  general.  One  expects  such  series 
to  represent  ^(x)  and  ce(x)  only  for  p  and  j3  that  characterize  non-percolating  states;  as 
yet  only  partial  results  are  available7'2*  for  their  radii  of  convergence.  We  note  that,  in 
general,  one  cannot  give  a  separate  physical  interpretation  in  terms  of  probability  densities, 
of  cc(x),  because,  unlike  gc(x),  it  need  not  be  positive  definite. 

The  expansion  procedure  of  Hill  can  be  carried  out  automatically  by  using  the  isomorphism23! 
between  percolation  and  the  one-state  limit  of  the  s-state  Potts  model.  Specifically,  the 
one-state  limit  of  a  continuum  Potts  model8  with  interparticle  potential 

Vij  =  <f>(xi ij)  +  v(*ij)[l  —  6ffit Tj]  (3.6) 

gives  a  correlated  continuum  percolation  model  with  interparticle  potential  <j>(x)  and 
separation- dependent  bond  probability  pj(x).  If  one  develops  Mayer  expansions  for  the 
thermodynamic  quantities  of  the  model  defined  by  (3.6),  and  applies  the  operator  —  |»=i 
to  them,  they  yield  the  basic  quantities  in  the  description  of  the  corresponding  percolation 
model.  This  procedure  provides  a  realization  of  the  general  percolation  process  described 
below  (3.5),  with  each  pair  of  particles  connected  with  a  separation-dependent  bond  prob¬ 
ability  given  by 

P6(*)  =  1  -  exp[-/3v(x)]  (3.7) 

with  u(x)  as  in  (3.6).  If  we  write  the  Ornstein-Zernike  equation  (2.1)  for  the  specific  case 
of  the  s-state  Potts  model,  takes  the  a  — *•  1  limit  and,  using  identities8,21 


h(x i,a,x2,0)  -»  -ge(x12) 

(3.8) 

c(xj,o,* 2,0)  -*  -cc(*  12) 

(3.9) 
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we  recover  the  equation  (3.1).  Here  a  and  fi  are  any  two  different  spin  states. 

When  we  restrict  the  continuum  Potts  model  to  a  lattice,  by  imposing  the  added 
restriction  that  particles  only  occupy  positions  whose  coordinates  are  integers,  the  result 
is  a  Potts  lattice  gas24  whose  one  6tate  limit  is  a  very  general  percolation  model.  Before 
doing  this,  we  add  to  the  potential  4>{x,j)  in  (13.6)  delta-function  potential  interaction 
that  prevents  overlap  of  two  lattice  gas  particles,  and  thus  of  two  sites  in  the  resulting 
percolation  model.  We  then  set  the  function  pfc(z)  equal  to  pb,  a  nearest-neighbor  bond 
probability  for  |x|  =  o,  and  equal  to  zero  for  x  ^  o.  The  lattice  site  and  bond  percolation 
models  are  given  by  special  cases  Pb  =  1,  and  p  =  1,  respectively. 

The  simplest  percolation  models  are  directly  related  to  thermal  lattice  models.  For 
example,  pure  bond  percolation  is  the  s  — *  1  limit  of  the  s-state  lattice  Pott6  model. 
Site  percolation  is  the  zero-temperature  limit  of  a  site-dilute  Ising  model.91  It  can  also 
be  realized  as  the  one-state  limit  of  a  Potts  model  containing  multi-Bite  interactions.31 
However,  the  former  mapping  seems  not  to  be  useful  in  this  context,  and  the  latter  has 
not  yet  been  exploited  in  this  context  because  of  its  complexity. 

We  now  use  these  mappings  to  construct  MSA-like  approximations  for  specific  perco¬ 
lation  models.  We  will  use  the  same  terminology  as  in  Section  2,  but  by  ho,  hlf  we  will 
mean  <?r(x)  evaluated  at  x  =  0  and  jz|  =  o,  respectively.  Similarly,  by  Co,  Ci,  we  denote 
ce(x)  for  x  =  0  and  |x|  =  a,  respectively.  Site  percolation  models  have  already  been  studied 
using  the  MSA6  approach.  In  this  case  it  is  natural  to  choose 


ho  =  0 


hi  =  1, 


(3.10) 


the  former  because  we  must  forbid  multiple  occupation  of  sites,  as  just  discussed,  and  the 
latter  because  neighboring  occupied  sites  are  always  connected.  Substituting  (3.9)  into 
(2.19)  gives  the  critical  site  density  for  percolation 


zQo  ~  1 

~  ”737" 


(3.11) 
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This  is  found  numerically  to  be  a  good  approximation  in  general.6 

For  hypercubic  lattices,  the  MSA  just  described  gives  pc  to  high  accuracy  for  d  >  4. 
However,  in  three  dimensions,  the  scheme  gives  pc  =  .341,  where  series  analysis  gives 
.312  ±  .003.8°  This  discrepancy  indicates  that  an  optimal  MSA  for  site  percolation  has 
yet  to  be  found.  We  note  that  the  effective  medium  approximation  (EMA)18  for  the 
conductivity  of  a  site-disordered  resistor  lattice  has  similar  difficulties.  The  EMA  is  very 
similar  in  spirit  to  the  hybrid  approximations  to  be  discussed  next.  Although  this  scheme 
gives  excellent  approximations  in  two  dimensions,  it  also  predicts14  a  three  dimensional 
site  percolation  threshold  which  is  too  low  by  10%. 

As  already  discussed,  continuum  percolation  has  many  similarities  to  site  percolation 
as  well  as  bond  percolation.  Thus  the  difficulty  just  discussed  may  also  account  for  the 
need  to  add  correction  terms8'5  to  the  naive  MSA  for  random  continuum  percolation  in 
order  to  recover  a  good  estimate  of  the  critical  point  in  this  model.  These  matters  are 
under  investigation. 

In  the  case  of  bond  percolation,  we  keep  the  first  equation  of  (3.10),  but  must  modify 
the  second.  One  possibility  is  to  follow  the  intuitive  notion  that  the  direct  connectedness 
function  c(*i2)  should  be  the  probability  density  associated  with  having  a  direct  bond 
between  the  sites  at  *i  and  xj 

ci  =  pb  (3.12) 

Using  the  Potts-model  correspondence,  and  keeping  only  terms  to  first  order  in  p  also  gives 
this  approximation.  This  form  of  the  the  MSA  is  appropriate  in  conditions  of  low  bond- 
density;  in  terms  of  the  related  Potts  model,  pb  is  also  a  ’’high-temperature”  expansion 
variable  (see  e.g.  3.7).  Also,  we  set  p  =  1  because  all  sites  are  occupied  in  bond  percolation. 
This  gives  for  the  bond  percolation  threshold 

(Pb)erit  =  &0  (3.13) 

Sahimi  et.al.38  noted  from  numerical  comparison  that  this  relation  provided  an  extremely 
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good  approximation  for  all  regular  lattices,  and  in  all  dimensions  d  >  3.  It  is  very  satisfying 
that  an  intuitively  reasonable  form  of  the  MSA  gives  just  this  result.  Table  1  shows 
the  quality  of  the  estimate  (3.13)  for  a  variety  of  three-dimensional  lattices,  as  well  as 
hypercubic  lattices  in  higher  dimensions. 

We  could  have  instead  made  the  approximation 

Ci  =  K  (3.14) 

where  K  is  defined  by  the  correspondence  between  the  Potts  model  coupling  constant,  and 
the  bond  probability  in  the  percolation  model  which  is  its  one-state  limit: 

pb  =  1  -  t~K  (3.15) 

In  other  words,  if  we  very  formally  attempt  to  treat  bond  percolation  as  if  it  were  a 
thermal  model  (remembering  the  Potts-model  correspondence),  we  would  be  led  to  try  the 
high-temperature  approximation  (3.14).  This  gives  the  critical  condition 

(Pb)cnt  =  1  -  exp(-Go)  ~  .223  (3.16) 

This  should  be  compared  with  (3.13),  which  gives  a  value  of  .252,  and  the  existing  series 
analysis  results,  which  give  .249  ±  .0002.  Thus  (3.16)  is  a  reasonable  approximation,  but  it 
lacks  the  remarkable  accuracy  of  (3.13).  A  similar  result  is  found5  in  studies  of  continuum 
percolation,  where  a  form  for  c(z)  must  be  assumed  over  the  entire  range  of  values  for 
which  v(z)  is  nonzero.  In  that  case  also,  the  low-density  ansatz  (3.12)  is  found  to  give  a 
prediction  for  the  threshold  which  is  numerically  superior  to  the  high-temperature  ansatz 
(3.14). 

We  make  the  observation  that  approximations  for  the  short-range  values  of  c(x)  would 
be  better  motivated  if  in  fact  that  quantity  were  a  probability  density;  in  fact,  it  seems  never 
to  be  positive  definite.  This  can  be  easily  checked  for  the  problems  studied  here  because 
c(z)  takes  only  two  values,  Co  and  e%.  For  both  pure  site  and  pure  bond  percolation,  it 
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is  always  found  that  for  any  lattice  studied  that  pc\  >  0  and  pco  <  0.  If  this  were  not 
true,  then  critical  percolation,  in  the  MSA,  would  be  directly  isomorphic  to  a  random  walk 
model26’27  defined  by  transition  probabilities  pc o,  pc\.  However,  since  having  a  nonzero 
value  for  Co  simply  rescales  the  time  coordinate  describing  the  progress  of  a  random  walk, 
one  can  always  rescale  the  other  nonzero  values  of  c,  by  (1  —  pc0)  and  get  a  physically 
realizable  random  walk.  In  fact,  we  can  rewrite  the  basic  equation  (2.15)  in  a  way  that  is 
applicable  to  MSA-like  approximations  in  which  c(x)  is  not  assumed  to  vanish  for  sr  ">  1: 


=  phi 


(3.17) 


(1  -pc0) 

The  existence  of  this  formal  equivalence  is  seen  to  be  a  general  fact  about  all  mean  spherical 
approximations,  even  those  for  anisotropic  or  directed  models,  as  we  discuss  in  Section  5. 
However,  in  general,  the  coefficients  c<  oscillate  in  sign,  thus  higher-order  approximations 
do  not  give  realizable  walks.  We  note  that  the  normalization  condition  for  the  transition 
rates  in  such  a  random  walk  is  just  the  criticality  condition  for  the  model  being  studied, 
(see  e.g.2.14). 

Since  the  MSA  gives  a  good  approximation  to  the  threshold  for  both  pure  site  and 
pure  bond  percolation,  it  is  natural  to  use  it  to  study  the  general  site-bond  percolation 
model,  in  which  a  cluster  is  defined  to  be  a  group  of  occupied  sites  connected  by  occupied 
bonds.  As  before,  there  are  several  natural  approximations  that  one  can  use  to  close  the 
Ornstein-Zernike  equation.  Note  that  (3.11)  is,  a  priori,  just  as  reasonable  an  assumption 
in  the  general  site- bond  problem  as  in  the  pure  bond  problem.  Using  it  in  the  general 
problem,  however,  gives  a  percolation  locus  in  the  (p,pb)  plane  defined  by 


PPb  =  Go 


(3.18) 


which,  e.g.,  in  the  case  of  pure  site  percolation,  is  immediately  seen  to  be  a  very  poor 
approximation. 

Thus,  we  instead  approximate  hi  by  enumerating  the  smallest  graphs  that  contribute 
to  it,  i.e.,  the  smallest  bond  sets  that  join  two  sites  that  are  nearest  neighbors.  This  is 
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equivalent  to  using  the  Potts  lattice-gas  mapping  described  below  (3.8)  and  calculating 
the  Mayer  expansion  of  that  model.  On  the  simple  cubic  lattic,  considering  just  the  two 
graphs  of  Figure  1  gives 


hi(p,Pb)  =  Pb  +  (1  -  Pb)  * 


1  -  (1  -p2pI)' 


(3.19) 


~  Pb  +  (1  -  Pb)Gp2Pb  +  0{pA) 

This  approximation  for  hi  is  exact  for  site  percolation,  and  gives,  for  bond  percolation,  the 
critical  value  ( hi)Crtt  =  <258,  whereas  the  approximation  based  on  (3.11)  gives  .252  for  the 
same  quantity.  The  percolation  locus  in  the  (p,Pb)  plane  as  given  by  the  approximation 
(3.16)  is  shown  in  Figure  2.  This  locus  has  been  obtained  by  simulation  in  both  two 
dimensions28  and  three  dimensions.29  The  critical  locus  given  by  substituting  (3.19)  in 
(2.19)  is  found  by  calculating  the  critical  site  density  p ,  for  a  specific  value  of  pj,,  by  using 
the  Newton-Raphson  method.  This  approximation  is  already  of  high  quality  and  can  easily 
be  improved  by  adding  terms,  except  for  the  part  of  the  phase  plane  near  the  pure  6ite 
percolation  limit;  we  discuss  this  problem  further  in  Section  5.  It  iB  worth  noting  that 
many  different  schemes  are  available  for  estimating  the  quantity  hi  in  both  thermal  and 
percolation  models.  In  particular,  the  method  of  Kikuchi31  involves  assuming  a  functional 
form  for  the  free  energy  which  contains  as  parameters  the  values  of  h(x)  for  small  separation 
x,  and  mimimizing  this  functional  to  determine  these  quantities.  We  discuss  this  class  of 
approximations  in  the  next  section. 

4.  Hybrids  of  the  Kikuchi  Cluster  Approximation  and  MSA 
In  Section  3,  we  showed  that  one  can  obtain  good  estimates  for  phase  transition  loci 
from  MSA-like  approximation  schemes  if  reasonable  estimates  are  available  for  values  of 
the  correlation  functions  at  short-range.  In  this  section,  we  explore  the  possibility  of 
using  the  Kikuchi  cluster  variational  method  (CVM)31  to  determine  these.  In  order  to 
complement  our  discussion  in  Section  3  of  the  basis  of  an  MSA  approach  to  percolation 
models,  we  also  obtain  directly  a  Kikuchi  CVM  for  bond  percolation,  by  using  the  Fortuin- 
Kastelyn  mapping23  between  percolation  and  the  Potts  model.  We  compare  the  direct 
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estimates  of  the  critical  bond  probability  given  by  this  method  with  hybrid  estimates  given 
by  using  the  structural  information  from  this  method  as  input  to  MSA-like  approximations. 
The  hybrid  estimates  are  found  to  be  superior.  The  cluster  variational  method(CVM), 
first  systematically  developed  by  Kikuchi81,  involves  several  steps.  We  outline  these  here, 
restricting  our  discussion  to  the  Ising-like  spin  models  for  which  the  theory  was  originally 
developed: 

1.  Because  Ising  variables  take  discrete  values,  the  various  small  clusters  of  contiguous 
spins  (pair6,  triplets,  etc.)  can  take  on  only  a  finite  number  of  possible  configurations. 
After  specifying  a  set  of  small  dusters  to  serve  as  a  basis  set,  one  chooses  as  working 
variables  the  probabilities  of  occurrence  of  each  possible  configuration  of  these  clusters. 
If  the  basis  set  consists  of  only  one  cluster,  a  nearest-neighbor  pair,  the  corresponding 
occurrence  probabilities  are  just  the  values  of  the  spin-spin  correlation  function  at  nearest- 
neighbor  separation.  We  note  that  these  variables  are  exactly  the  quantities  needed  in 
the  MSA-like  approximations  discussed  in  this  paper.  Direct  use  of  these  values  yields  the 
well-known  Bethe  approximation.24 

2.  In  terms  of  the  working  variables,  one  writes  a  consistent  approximation  for  the  free 
energy  of  the  system.  Requiring  that  this  expression  be  minimized  with  respect  to  the 
working  variables  then  gives  a  set  of  constraint  equations  to  determine  the  values  of  the 
working  variables  as  functions  of  the  system  parameters  (temperaure,  magnetic  field,  etc.) 

3.  One  can  then  determine  an  approximation  for  the  critical  point,  by  requiring  that  the 
symmetry-breaking  variables  display  power-law  behavior  as  the  singularity  is  approached. 

Instead  of  following  this  last  step  to  determine  the  location  of  the  critical  point  one  can 
instead  use  the  following  hybrid  method:  use  the  functional  expressions  for  values  of  the 
two-point  correlation  function  at  small  separation  as  input  to  the  MSA-like  approximations 
discussed  in  Section  3.  Sperifically,  the  CVM  will  give  expressions  for  the  quantities  ho 
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and  hi  in  the  equation 


(4.1) 


phi  zGq  —  1 
1  +  pho  zGo 

as  given  by  the  MSA.  We  remark  that  it  is  not  clear  a  priori  that  this  method  of  determining 
the  critical  point  will  be  successful,  as  no  single  consistent  expansion  scheme  has  been 
employed.  For  example,  the  numbers  used  to  evaluate  the  RHS  of  (4.1)  are  obtained  by 
setting  c(x)  to  zero  for  |x|  >  a;  this  type  of  low-density  approximation  is  not  a  priori 
consistent  with  use  of  an  extremely  accurate  value  of  h\. 

The  algebra  involved  in  realizing  the  program  outlined  above  has  been  detailed  in 
the  beautiful  paper  of  Kikuchisl°  and  will  not  be  repeated  here.  We  follow  the  notation 
used  in  that  paper  and  merely  give  the  results  of  our  calculations.  The  variable  hi  used 
in  Section  2  to  describe  nearest-neighbor  values  of  the  Ising  model  correlation  function  is 
related  to  Kikuchi’s  variable  yi  by  hi  =  4j/i  —  1.  If  we  use  the  lowest-order  CVM,  in  which 
the  only  cluster  in  the  basis  set  is  a  nearest-neighbor  pair,  the  result  is 

H2  -  H~2 
hl  ~  6  -  H2  -  H~2 

with  H  =  cxp[K ising]-  This  result  was  calculated  for  a  two-dimensional  Ising  model,  but 
at  this  low  level  of  approximation,  it  is  entirely  consistent  to  use  the  relation 

^Rising  =  Rising  (4-3) 

noting  that  the  Bethe  approximation  per  se  depends  only  on  the  combination  zK,  where 
z  is  the  coordination  number.  Substituting  (4.2)  into  (4.1)  gives  a  critical  point  located 
at  Kerit  ~  -779,  as  compared  with  the  value  Kcrn  =  .918  given  by  series  analysis.  This 
is  slightly  better  than  the  value  Keru  =  1.099  given  by  the  direct  Kikuchi  method  at  the 
same  level  of  approximation.  An  improved  treatment  which  makes  explicit  use  of  the  4D 
nature  of  the  lattice3 lB  gives 

h  +*-1 

1  +  6<f>  +  1 


20 


where  the  auxilliary  variable  4>  is  defined  implicitly  by 

B2  =  exp  [^Klg] 

_  1 
"  <t> 

Substituting  (4.4)  into  (4.1)  gives  Kcra  =  .926,  which  is  an  extremely  good  approximation! 

In  order  to  apply  the  same  approximation  scheme  to  percolation  models,  we  first 
develop  the  Kikuchi  cluster  approximation  for  bond  percolation.  Kikuchi31  applied  his 
method  to  site  percolation  by  noting  that  it  is  equivalent  to  the  zero-temperature  Ising 
model.  He  treated  bond  percolation  as  site  percolation  on  the  corresponding  alternate 
lattice.  For  many  common  lattices,  e.g,  the  simple  cubic  latice,  this  requires  an  enlarged 
primitive  cell  and,  presumably,  requires  including  larger  clusters  in  the  basis  set  to  give 
results  of  comparable  accuracy  to  that  obtained  for  site  percolation.  We  proceed  instead 
by  calculating  KCTn(s)  using  the  CVM  for  a  dilute  s-state  Potts  model,  then  taking  the 
one-state  limit  as  described  in  Section  3.  Using  the  Bethe  approximation  for  the  two- 
dimensional  square  lattice  gives  for  the  bond  percolation  threshold  =  .4226,  as  compared 
with  the  exact  result  p&  =  .5.  Here  we  used  the  correspondence  (3.15)  between  the  Potts- 
model  coupling  and  bond  probability. 

We  now  estimate  the  three-  dimensional  bond  percolation  transition  by  using  a  pro¬ 
cedure  parallel  to  that  used  above  for  the  Bethe  approximation  to  the  Ising  model.  First 
we  use  the  direct  Kikuchi  procedure  just  described  which  is  based  on  the  Potts-model 
mapping.  Using  the  scaling  6 =  4 as  above  gives  pb  =  .231  for  the  transition 
point  in  a  simple  cubic  lattice,  as  compared  to  the  value  .252  given  by  series  methods. 
We  now  instead  use  the  functional  form  for  the  nearest-neighbor  connectedness  function, 
as  given  by  the  Kikuchi  method  just  described,  as  input  to  the  MSA  defined  by  (4.1). 
The  Bethe  approximation  for  the  Potts-model  variable  yi ,  which  is  the  probability  that  a 
nearest-neighbor  pair  are  in  different  spin  states,  is  given  by 

Vi  =  l./l*  +  »(»  ~  l)exp(-2 K)\  (4.6) 


3  <f>  -t-  1 
^  +  3 
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We  obtain  the  nearest-neighbor  connectedness  function  h\  for  bond  percolation  by  using 
the  correspondence  (3.8)  and  the  relation  hi  =  4yi  -  1.  Using  the  result  in  the  MSA  as 
before  gives  a  critical  bond  probability  pj,  =  .242,  a  substantially  more  accurate  result. 

As  a  more  involved  illustration  of  this  hybrid  Kikuchi-MSA  method,  we  consider  the 
Wannier  approximation,  in  which  the  nearest-neighbor  pair  and  elementary  plaquette,  or 
square,  are  taken  to  be  basic  clusters.  It  seems  quite  difficult  to  use  the  direct  Kikuchi 
method  described  above,  in  which  one  first  applies  the  cluster  variational  method  to  the 
s-state  Potts  model,  then  takes  the  one-state  limit,  as  a  means  for  locating  the  percolation 
critical  point.  Indeed,  if  the  spins  in  the  basic  clusters  are  allowed  to  be  in  any  of  a 
states,  with  one  of  the  states  distinguished  as  the  symmetry-breaking  state,  there  are  20 
different  configurations  of  these  clusters  (see  Figure  3).  Locating  the  critical  point  then 
involves  finding  the  determinant  of  a  matrix  of  rank  20.  We  note  in  passing  that  Kikuchi’s 
method81  of  realizing  bond  percolation  as  site  percolation  on  an  alternate  lattice  is  easily 
seen  to  lead  to  an  equally  large  variable  set.  The  hybrid  MSA  technique  developed  in  this 
section  is  readily  applied  to  this  model,  however.  Since  the  one-stt.te  limit  of  the  Kikuchi 
method  is  of  some  interest  in  its  own  right,  we  describe  it  in  Appendix  A.  Here  we  sketch 
the  procedure  and  discuss  the  results  of  using  it.  After  setting  the  number  of  states  a 
equal  to  one,  the  values  of  the  configuration  variables  in  the  symmetric  state  are  quickly 
obtained.  We  derive  an  expression  for  the  pair  configuration  variable  yu  as  a  function  of 
bond  density,  along  with  the  corresponding  form  for  hi.  Substituting  this  into  the  MSA 
equation  (4.1)  gives  a  percolation  critical  point  at  bond  density  pe  =  .249.  This  is  in 
remarkable  agreement  with  the  best  series  estimate,32  pc  =  .249  ±  .0002! 

5.  MSA-Like  Approximations  for  Directed  Percolation 

In  this  section  we  apply  the  class  of  MSA-like  approximations  discussed  in  this  paper 
to  directed  site  and  bond  percolation  models. 

In  order  to  gain  some  perspective  on  the  strengths  and  limitations  of  the  methods 
discussed  in  this  paper,  we  use  them  to  calculate  the  two-point  connectedness  function, 
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and  the  location  of  the  critical  point,  in  directed  percolation  models.  We  can  define  these 
in  general  as  follows:  a  distinguished  direction  is  chosen  in  the  space  occupied  by  the  lattice 
being  studied.  This  direction  may  or  may  not  coincide  with  one  of  the  principle  axes  of  the 
lattice.  When  a  vector  in  this  preferred  direction  is  projected  onto  the  bonds  of  the  lattice 
it  induces  in  them  an  allowed  direction  of  passage.  In  the  convention  adapted  here,  bonds 
which  are  orthogonal  to  the  preferred  vector  remain  non-oriented  and  thus  allow  two-way 
passage  or  connection.  Models  containing  a  class  of  such  bonds  will  then  be  called  partly 
directed  percolation  models.  This  construction  is  motivated  by  one  of  the  basic  classes 
of  applications  for  directed  percolation  models:  transport  through  random  or  two-phase 
materials  under  the  influence  of  a  uniform  gradient  or  bias  field.  A  given  lattice  may  then 
yield  a  number  of  different  directed  or  partly  directed  models  depending  on  the  preferred 
direction  chosen.  In  terms  of  MSA-like  methods,  these  models  differ  from  isotropic  models 
in  one  basic  way:  random  walks  on  the  corresponding  lattices  either  allow  only  a  very 
restricted  class  of  closed  paths,  or  allow  no  such  paths.  The  importance  of  this  fact  will 
be  explored  further  when  we  discuss  the  results  of  this  section  in  general  terms.  Here  we 
focus  on  developing  specific  MSA-like  approximations. 

We  may  directly  adapt  the  methods  already  developed  to  treat  directed  percolation  as 
regular  percolation  on  a  lattice  with  peculiar  conductivity  properties.  Specifically,  we  make 
equations  (2.9)  and  (4.1)  the  basis  of  our  treatment.  The  latter  will  be  used  unchanged, 
while  for  the  SC  directed  lattice  (with  the  (1,1,1)  vector  the  preferred  direction),  the 
equation  (2.9)  becomes 


d*k 


(5.1) 


x  n  f _ _ 

*,o  +  P  (*  “(2 7r)s  y  1  —  pc0  —  c(xj) exp(fc  •  Xj) 

=  <?(*) 

Here  the  sum  is  over  all  sites,  at  positions  x* ,  for  which  Cj  /  0.  In  general,  directed 
percolation  models  will  have  two  correlation  lengths,  both  of  which  become  infinite  at  the 
critical  point.*2  We  focus  here  on  the  value  predicted  by  MSA  for  the  critical  point.  We 
calculated  n**>order  approximations  as  follows: 
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1.  We  assume  that  the  n  values  of  c(x)  corresponding  to  nearest  neighbor,  next-nearest 
neighbor,  etc.,  separation  are  nonzero;  for  larger  separations  x  we  assume  c(x)  =  0.  The 
nonzero  values  c,  are  related  by  the  random- walk  representation  (5.1)  to  nonzero  values  of 
step  fugacities  by 

Wi  =  c,  (5.2) 

Here  «>i,  tnj,  etc.  are  fugacities  for  steps  to  nearest-neighbor,  next-nearest  neighbor,  etc., 
sites.  We  have  used  the  fact  that  Co  =  0  for  this  model;  this  can  be  seen  by  setting  *  =  0 
in  (5.1). 

2.  One  then  solves  n  equations  of  form  (4.1)  for  the  step  fugacities  w<  as  functions  of  the 
(site  or  bond)  density.  Here  we  use  explicit,  exact  formulas  for  both  the  connectedness 
function  hi  and  the  random- walk  generating  function  zG{.  These  are  readily  found  because 
the  directed  random  paths  between  any  pair  of  points  which  contribute  to  these  functions 
are  self-avoiding  walks,  i.e.,  they  lack  loops.  Also,  the  number  of  such  walks  i6  small. 

3.  Substituting  the  exact  fugacities  Wi(p)  into  the  criticality  condition  =  1  then 

gives  a  polynomial  equation  whose  smallest  positive  root  is  the  critical  density. 

The  critical  densities  given  by  successive  approximations  of  this  type  are  listed  in 
Table  2.  These  numbers  seem  to  converge;  however,  for  both  site  and  bond  problems,  the 
resulting  critical  density  is  lower  than  the  simulation  value  by  about  8%. 

6.  Limitations  of  the  MSA  and  Directions 
for  Further  Research 

In  this  section,  we  analyze  possible  reasons  for  the  failure  of  MSA-like  methods  to 
yield  highly  accurate  percolation  thresholds  for  some  systems. 

Why  do  MSA-like  methods  give  substantially  better  threshold  values  for  bond  per¬ 
colation  models  than  for  site  percolation  models?  Of  course,  the  extension  of  MSA-like 
methods  in  the  former  case  is  better  motivated  than  in  the  latter  case;  this  was  the  pur¬ 
pose  of  our  development  of  the  Potts-model  formalism  in  Section  3.  But  we  need  a  more 
basic  understanding  to  extend  these  methods  further.  Here  we  will  explore  two  possible 
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elements  in  such  an  understanding. 

The  dominant  singularity  in  the  Mayer  expansion  of  mean  cluster  size,  and  other  phys¬ 
ical  quantities  that  describe  bond  percolation  models,  is  the  physical  percolation  threshold. 
This  is  not  true  in  general  for  site  percolation  models;  their  dominant  singularity  tends  to 
be  located  at  a  negative  retd  value  of  density.  Dominance  by  unphysical  singularities  has 
also  been  found  in  the  series  expansions  of  directed  site  models.32 

We  have  no  general  argument  that  MSA-like  approximations  will  yield  a  real,  positive 
density  as  the  dominant  singularity.  However,  this  is  found  to  be  the  case  with  all  the 
approximations  studied  in  this  paper.  The  dominant  singularity  in  the  anti-ferromagnetic 
lattice  Potts  model  lies  on  the  negative  real  axis.  In  the  MSA,  the  Potts  model  for  neg¬ 
ative  density  is  mapped  onto  bond  percolation  at  positive  density;  thus,  in  the  MSA,  the 
dominant  singularity  of  bond  percolation  occurs  at  a  positive,  physical  density.  We  ob¬ 
serve  the  same  fact  in  the  MSA-like  approximations  studied  in  Section  5.  If  one  plots  the 
singularities  of  the  mean  cluster  size,  which  are  just  the  zeroes  of  the  polynomial  equation 

pc(k  =  0,p)  =  1  (6.1) 

one  always  finds  the  dominant  singularity  at  a  physical  density.  The  unphysical  singu¬ 
larities  associated  with  site  percolation  also  cause  difficulties  in  applying  other  standard 
methods  for  studying  phase  transitions.  Dealing  with  such  structures  via  approximations 
of  MSA  type  is  thus  an  open  problem. 

The  slow  convergence  exhibited  by  the  MSA  for  3D  directed  problems  can  be  under¬ 
stood  in  two  complementary  ways.  We  briefly  describe  both  of  them.  In  the  direct  form  of 
the  MSA  described  in  Section  5  for  directed  percolation  models,  the  connectedness  func¬ 
tion  is  represented  as  a  generating  function  for  directed  random  walks.  The  MSA,  roughly 
speaking,  uses  the  balance  between  random  walks  that  return  to  their  starting  point  and 
those  that  do  not  to  capture  the  balance  between  the  short  range  and  long-range  behavior 
of  ge( x)  at  criticality.  In  directed  models  in  which  all  walks  lack  recurrences,  this  balance 
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is  missing  because  only  short  walks  can  contribute  to  the  approximations  for  gc(z)  at  small 
separation. 

Equivalently,  one  can  reformulate  the  problem  of  calculating  the  percolation  threshold 
in  terms  of  random  walks  with  recurrence,  but  in  (e£  —  1)  dimensions.  If  a  directed  3D 
percolation  cluster  is  projected  onto  the  plane  perpendicular  to  the  preferred  direction,  each 
site  can  be  identified  with  a  two-dimensional  vector  r_L ,  its  position  in  the  perpendicular 
plane  with  respect  to  the  origin  of  that  cluster.  Then  define*2 


5(f±)  =  J  dr||$c(r,|,rj.  =  fi_) 

=  J  dk±  cxp[-ik±  •  r±]ge(kj_,k\\  =  0) 


(6.2) 


to  be  the  expected  number  of  sites  contained  in  the  cluster  and  located  at  lateral  positions 
f _L-  Proceeding  as  in  Section  3  then  shows  that  is  given  by  the  generating  function 

of  a  random  walk  process  in  ( d  —  1)  dimensions.  This  process  occurs,  in  general,  on  a 
directed  lattice  (for  the  3D  SC  lattice  the  corresponding  process  occurs  on  the  2D  directed 
cyclic  trianglular  lattice),  but  involves  random  walks  with  nonzero  probability  of  return  to 
the  origin.  However,  this  formulation  shows  that  the  MSA  describes  a  directed  3D  process 
in  terms  of  a  2D  process.  Since  the  MSA  is  in  general  inapplicable  in  two  dimensions,  this 
gives  another  view  of  the  failure  of  the  MSA  in  this  case.  From  this  analysis,  however,  one 
expects  the  MSA  to  give  accurate  critical  densities  for  higher-dimensional  (>  3)  directed 
percolation. 

7.  Conclusions 

The  MSA  approach  to  site  percolation  has  been  successfully  extended  to  both  bond 
percolation  and  general  site-bond  percolation.  The  bond  percolation  threshold  given  by 
this  method  is  found  to  coincide  with  an  analytic  estimate  already  shown  to  be  of  high 
accuracy.  For  the  general  site-bond  percolation  model  on  the  SC  lattice,  the  percolation 
locus  calculated  from  this  approximation  agrees  quite  well  with  that  given  by  simulation. 
Better  agreement  will  require  a  general,  reliable  method  for  treating  site  percolation  mod- 
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els,  possibly  using  the  corresponding  Potts-model  mapping.22 

A  class  of  approximations  of  MSA  type  have  been  applied  to  directed  site  und  bond 
percolation  models.  These  approximations  give  moderate  accuracy,  but  will  require  ba¬ 
sic  improvements  to  give  highly  accurate  predictions.  Some  reasons  for  this  have  been 
identified. 

A  major  advantage  of  this  approach  to  percolation  is  that  the  resulting  integral  equa¬ 
tions  can  be  solved  analytically  to  give  the  connectedness  function  £e(xi2)*  Summing  this 
function  over  all  possible  separation  vectors  £12  then  gives  the  mean  cluster  size.  The 
results  described  in  this  paper  can  be  generalized  in  a  number  of  directions  without  giving 
up  this  advantage.  For  example,  the  bond  probability,  which  in  this  paper  was  taken  to 
be  non-zero  only  for  particles  with  nearest-neighbor  separation,  can  be  taken  to  have  cer¬ 
tain,  non-trivial,  long-range  forms  while  still  allowing  exact  solution  for  the  connectedness 
function. 

It  would  be  valuable  to  have  efficient  computational  procedures  for  the  accurate  de¬ 
termination  of  gc(xi  2)  in  a  general  correlated  percolation  problem.  The  hybrid  procedure 
discussed  in  Section  4,  in  which  the  Kikuchi  method  is  used  to  calculate  the  short-range 
values  of  ge(*12),  and  these  are  then  used  in  the  MSA,  shows  great  promise  in  prelim¬ 
inary  studies  reported  here.  It  would  be  useful  to  find  a  direct  Kikuchi  approximation 
for  bond  percolation,  so  as  to  eliminate  the  added  algebraic  complexity  introduced  by  the 
Potts-model  map,  if  this  is  possible. 

Also,  in  some  models  it  may  be  necessary  to  use  clusters  substantially  larger  than 
those  tractable  by  analytic  means.  For  these  models,  an  analog  to  the  numerical  methods 
used  in  the  phenomenological  renormalization  group  would  be  valuable. 

These  matters  are  presently  under  study. 
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In  this  appendix,  we  sketch  the  Wannier,  or  four-site,  approximation  for  bond  perco¬ 
lation  on  the  simple  cubic  (SC)  lattice. 

In  the  Wannier,  or  four-site,  approximation,  the  bond  and  square  composed  of  near¬ 
est  neighbor  sites  are  chosen  as  the  set  of  small  clusters  used  as  a  basis  for  building  up 
correlations.  For  a  general,  s-state  Potts  model,  there  are  20  different  configurations  of  the 
spins  in  these  clusters.  These  are  shown  in  Figure  3.  The  variables  giving  the  probability 
of  occurrence  of  each  site,  bond  and  square  configuration  are  denoted  by  *,•,  ,  and  Zijki, 

respectively,  where  the  subscripts  give  the  values  of  the  site  variables  involved.  In  terms 
of  these  variables,  the  Kikuchi  method  gives,  for  the  free  energy  of  the  system: 

9  E  Ivy ln  V*i  ~  Vii\  ~  7  E  ln  ~  3  S  t ZW ln  Ziikl  ~  «J 

*  j  «  ijkl 


c 

-X> 

E  yi>  ~ Xi 
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E  Zi’ki  ~  w 

.  t 

t 

j 

*j 

.  kl 

Here  z  is  the  coordination  number  and  K  is  the  Potts  coupling  constant.  The  first  term  in 
this  expression  is  the  energy  per  spin,  the  term  ir  urly  brackets  gives  the  entropy  per  spin, 
and  the  last  three  terms  incorporate  constraints  due  to  the  normalization  of  variables. 
The  expression  (A.l)  is  minimized  under  the  conditions 


}  (a.i) 


7  ln  Xi  —  C  +  D{  =  0 


(A.  2) 


\zK{*  -  6ij)  -  91n ytj  -  Di  +  Eij  -  0 


(A- 3) 


E 

4c  P. 


3  ln  Zijki  —  Eij 


=  0 


(AA) 


where  the  sum  in  the  equation  (A.4)  is  over  the  four  cyclic  permutations  of  {ijkl).  It  is  not 
difficult  to  solve  these  equations  in  the  symmetric  phase.  In  the  limit  s  — ►  1,  the  variables 
with  all  indices  equal  take  the  value  unity;  in  this  limit  they  give  the  correlations  of  a  non¬ 
interacting  spin  model.  The  one-state  limit  of  variables  whose  indices  are  not  all  equal 
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is  also  well-defined,  however,  and  gives  the  correlation  functions  of  bond  percolation.7’6 
For  example,  the  one-state  limit  of  yu  gives  the  nearest-neighbor  value  of  the  blocking 
function;  this  is  the  lattice  analog  of  (3.8).  This  quantity  is  given  by 


(A.  5) 


Vl2  =*V(1  -Pfc) 

=  [3®2  -  3x8  +  x4]8 

where  the  auxiliary  variable  x  =  exp[j^£12]  with  E12  the  Lagrange  multiplier  in  (A.l). 
Equation  (A.5)  determines  the  nearest-neighbor  corelation  function  hi  =  1  —  yn  as  a 
function  of  the  bond  density  p&.  Substituting  this  into  (2.19)  gives,  for  the  bond  percolation 
threshold  in  the  SC  lattice,  pb  =  .246,  which  is  in  excellent  agreement  with  the  series 
estimate  pb  =  .249  ±  .0002 .2 5 
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FIGURE  CAPTIONS 

1.  The  lowest-order  graphs  in  the  density  expansion  of  hi,  the  nearest-neighbor  connect¬ 
edness  function,  as  given  by  (3.19).  This  quantity  is  required  in  the  MSA  for  a  general 
site-bond  percolation  model. 

2.  The  percolation  locus  for  random  site-bond  percolation  on  the  3D  simple  cubic  (SC) 
lattice,  as  given  both  by  simulation  (solid  line),  and  by  the  MSA  of  Section  3  (dotted  line). 

3.  Configurations  whose  probabilities  of  occurrence  form  the  working  variables  for  the 
Wannier  or  square  approximation  to  the  properties  of  the  s-state  Potts  model.  By  conven¬ 
tion,  ’1’  here  denotes  the  symmetry-breaking  state;  ’2’,  ’3’,  etc.,  denote  any  other  distinct 
states.  The  a  =  1  limit  of  these  variables  give  the  two  and  four-point  correlation  functions 
of  bond  percolation  at  small  separation. 
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Table  Captions 

1.  Approximations  to  the  bond  percolation  threshold  given  by  the  low  density  MSA  dis¬ 
cussed  in  Section  3.  The  sources  for  the  numerical  estimates  are:  aAdler  et.  al.  (1990),  h 
Sykes  et.  al.  (1976),  'Vyssotsky  et.  al.  (1961).  This  table  was  adapted  from  reference  25 
and  revised. 

2.  Approximations  to  the  3D  directed  site  and  bond  percolation  thresholds  given  by  the 
MSA  scheme  discussed  in  Section  4.  The  slow  convergence  is  believed  to  be  a  result  of  the 
two-dimensional  nature  of  the  problem  in  the  MSA. 
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Table  1.  Approximations  to  Bond  Percolation  Threshold  Given  by  (3.18) 


Lattice 

series  expansions 

Go 

simple  cubic 

0.2488  ±  0.0002° 

0.25273 

body-centered  cubic 

0.18025  ±  0.00015“ 

0.17415 

face-centered  cubic 

0.119  ±0.001* 

0.11206 

hexagonal  dose-packed 

0.124  ±  0.005c 

0.11206 

4D  simple  cubic 

0.16005  ±  0.00015° 

0.156 

5D  simple  cubic 

0.11819  ±0.00004“ 

0.115 
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Table  2.  Approximations  to  Directed  Percolation  Threshold  Given  by  Equation  (3.18) 


Order  of  Approximation  Bond  Percolation  Site  Percolation 


1 

0.333 

0.333 

2 

0.348 

0.395 

3 

0.352 

0.396 

4 

0.356 

0.396 

5 

0.362 

0.396 

(EXACT) 

0.384 

0.432 

3 


